Następny etap to wygenerowanie układu równań liniowych oraz jego rozwiązanie.
Oczywiście zależy nam na tym żeby wybrane algorytmy generacji oraz rozwiązywania układów równań były najszybsze możliwe.
W pierwszej kolejności generujemy macierz zawierającą całki z iloczynów dwuwymiarowych funkcji B-spline.
Możemy dla uproszczenia przyjąć, że wektory węzłów są równoodległe, oraz, że jednostką odległości wzdłuż danej osi układu współrzędnych jest odległość pomiędzy dwoma węzłami w wektorze węzłów.
Zauważmy najpierw, że
\( \int B^x_{i,p}B^y_{j,p}B^x_{k,p}B^y_{l,p}dxdy = \int B^x_{i,p}B^x_{k,p}dx \int B^y_{j,p}B^y_{l,p}dy \)
Innymi słowy dzięki temu, że nasze dwuwymiarowe B-spline'y powstają poprzez przemnożenie jednowymiarowych B-spline'ów, możemy rozdzielić nasze dwuwymiarowe całki na iloczyn dwóch jednowymiarowych całek, w których przemnażamy i całkujemy pary jednowymiarowych B-spline'ów.
Musimy więc wygenerować następującą macierz

Zauważmy, że macierz ta ma następującą strukturę. W wierszach zmieniają się pierwsze B-spline'y w obu całkach pojedynczych (całce po $x$ i całce po $y$). Wiersze zmieniają się tak, że najpierw mamy kolejno B-spline'y
\( B^x_{1,p},B^x_{2,p},...,B^x_{N_x,p} \) dla ustalonego \( B^y_{1,p} \). Następnie powtarzamy wszystkie \( B^x_{1,p},B^x_{2,p},...,B^x_{N_x,p} \) dla ustalonego \( B^y_{2,p} \), i tak dalej, występują bloki wszystkich kolejnych B-spline'ów względem \( x \) dla kolejnych ustalonych B-spline'ów \( B^y_{j,p} \) aż do ostatniego bloku, w którym występują kolejno \( B^x_{1,p},B^x_{2,p},...,B^x_{N_x,p} \) dla ustalonego ostatniego \( B^y_{N_y,p} \)
W kolumnach sytuacja jest podobna, dotyczy ona natomiast drugich B-spline'ów w każdej całce pojedynczej. Kolumny zmieniają się więc tak, że najpierw mamy kolejno B-spline'y \( B^y_{1,p},B^y_{2,p},...,B^y_{N_y,p} \) dla ustalonego \( B^x_{1,p} \). Następnie powtarzamy wszystkie \( B^y_{1,p},B^y_{2,p},...,B^y_{N_y,p} \) dla ustalonego \( B^x_{2,p} \), i tak dalej, występują bloki wszystkich kolejnych B-spline'ów względem \( y \) dla kolejnych ustalonych B-spline'ów \( B^x_{i,p} \) aż do ostatniego bloku, w którym występują kolejno \( B^y_{1,p},B^y_{2,p},...,B^y_{N_y,p} \) dla ustalonego ostatniego \( B^x_{N_x,p} \)
Macierz, która ma taką strukturę jest iloczynem Kroneckera dwóch mniejszych macierzy, jednej, w której występują same wyrazy z pierwszych całek pojedynczych (pary B-spline'ów względem \( x \)),
i drugiej, w której występują same wyrazy z drugich całek pojedynczych (pary B-spline'ów względem \( y \)). Matematyczny zapis takiego iloczynu Kroneckera to symbol \( \otimes \) pomiędzy dwoma wspomnianymi macierzami. Symbol "=" z przodu oznacza iż iloczyn Kroneckea tych macierzy daje naszą dużą macierz.
\( =\begin{bmatrix} \int{B^x_{1,p}B^x_{1,p}}dx & \int{B^x_{1,p}B^x_{2,p}}dx & \cdots & \int{B^x_{1,p}B^x_{N_x,p}}dx \\ \int{B^x_{2,p}B^x_{1,p}}dx & \int{B^x_{2,p}B^x_{2,p}}dx & \cdots & \int{B^x_{2,p}B^x_{N_x,p}}dx\\ \vdots & \vdots & \vdots & \vdots \\ \int{B^x_{N_x,p}B^x_{1,p}}dx & \int{B^x_{N_x,p}B^x_{2,p}}dx & \cdots & \int{B^x_{N_x,p}B^x_{N_x,p}}dx \\ \end{bmatrix} \otimes \\ \begin{bmatrix} \int{B^y_{1,p}B^y_{1,p}}dy & \int{B^y_{1,p}B^y_{1,p}}dy & \cdots & \int{B^y_{1,p}B^y_{N_y,p}}dy\\ \int{B^y_{1,p}B^y_{1,p}}dy & \int{B^y_{1,p}B^y_{1,p}}dy & \cdots & \int{B^y_{1,p}B^y_{N_y,p}}dy \\ \vdots & \vdots & \vdots & \vdots \\ \int{B^y_{N_y,p}B^y_{1,p}}dy & \int{B^y_{N_y,p}B^y_{1,p}}dy & \cdots & \int{B^y_{N_y,p}B^y_{N_y,p}}dy \\ \end{bmatrix} \)
Zapis \( \otimes \) oznacza, iż przemnażamy przez siebie odpowiednie wyrazy dwóch macierzy tak żeby dostać oryginalną macierz ( 1 ).

Musimy teraz obliczyć całki z iloczynów par jednowymiarowych B-spline'ów, przedstawionych na Rys. 1.
Zauważmy najpierw, że jeśli przemnożymy przez siebie dwa jednowymiarowe B-spline'y, których wykresy nie nachodzą na siebie (matematycy powiedzą "które nie mają wspólnego suportu"), na przykład \( N_{1;2} \) oraz \( B_{4;2} \) z Rys. 1, wówczas całka
\( \int_{[t_0,t_6]} B_{1;2}(x)B_{4;2}(x)dx = 0 \)
Dzieje się tak dlatego ponieważ na każdym przedziale przedstawionym na Rys. 1 jeden z B-spline'ów jest równy zeru, czyli
\( \int_{[t_0,t_6]} B_{1;2}(x)B_{4;2}(x)dx = \\ \int_{[t_0,t_1]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_1,t_2]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_2,t_3]} B_{1;2}(x)B_{4;2}(x)dx +\\ \int_{[t_3,t_4]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_4,t_5]} B_{1;2}(x)B_{4;2}(x)dx + \int_{[t_5,t_6]} B_{1;2}(x)B_{4;2}(x)dx =\\ \int_{[t_0,t_1]} B_{1;2}(x)*0dx + \int_{[t_1,t_2]} B_{1;2}(x)*0dx + \int_{[t_2,t_3]} B_{1;2}(x)*0dx + \\ \int_{[t_3,t_4]} 0*B_{4;2}(x)dx + \int_{[t_4,t_5]} 0*B_{4;2}(x)dx + \int_{[t_5,t_6]} 0*B_{4;2}(x)dx = \\ 0+0+0+0+0+0=0 \)
Następnie, możemy przyjąć jako jednostkę długości średnicę elementu, wówczas każdy przedział \( [t_i,t_{i+1}] \) będzie miał średnicę równą jeden.
Następnie musimy wprowadzić wzory na poszczególne segmenty B-spline'a (przedstawione na prawym panelu Rys. 1 ). Segmenty te oznaczamy \( B_1,B_2,B_3 \).
\( B_1 (x)=\frac{1}{2} x^2 \\ B_2 (x)=-x^2+x+\frac{1}{2} \\ B_3 (x)=\frac{1}{2}(1-x)^2 \)
Ponieważ jedynie całki z nachodzących na siebie B-spline'ów dają niezerowe wartości, nad każdym przedziałem \( [t_i,t_{i+1}] \) musimy tak naprawdę obliczyć następujące całki, występujące w macierzy elementowej, w których nachodzą na siebie poszczególne segmenty B-spline'ów. Ponieważ podane wzory na segmenty obowiązują nad przedziałem 0,1, a całka nie zmienia się podczas przesunięcia (podobnie jak objętość pod stołem nie zmieni się jeśli przesuniemy stół na inne miejsce na płaskiej podłodze), możemy dla uproszczenia policzyć nasze całki właśnie w przedziale 0,1.
\( \begin{bmatrix} \int_0^1 B_1(x)B_1(x)dx & \int_0^1 B_1(x)B_2(x)dx & \int_0^1 B_1(x)B_3(x)dx \\ \int_0^1 B_2(x)B_1(x)dx & \int_0^1 B_2(x)B_2(x)dx & \int_0^1 B_2(x)B_3(x)dx \\ \int_0^1 B_3(x)B_1(x)dx & \int_0^1 B_3(x)B_2(x)dx & \int_0^1 B_3(x)B_3(x)dx \\ \end{bmatrix} \)
Innymi słowy mamy tylko 9 możliwości nachodzenia na siebie poszczególnych segmentów.
Ponieważ segmenty \( B_1 \) i \( B_3 \) są symetryczne, dodatkowo mamy symetrię macierzy, czyli tak naprawdę musimy policzyć sześć całek
\( \begin{bmatrix} \int_0^1 B_1(x)B_1(x)dx & \int_0^1 B_1(x)B_2(x)dx & \int_0^1 B_1(x)B_3(x)dx \\ == & \int_0^1 B_2(x)B_2(x)dx & \int_0^1 B_2(x)B_3(x)dx \\ == & == & \int_0^1 B_3(x)B_3(x)dx \\ \end{bmatrix} \)
Z symetrii segmentów wynika również iż segment \( B_3 \) ma kształt odwróconego segmentu \( B_1 \),
i dlatego \( \int_0^1 B_1(x)B_1(x)dx =\int_0^1 B_3(x)B_3(x)dx \), innymi słowy odwrócenie funkcji nie zmienia wartości całki (podobnie jak przekręcenie stołu o 180 stopni nie zmienia objętości znajdującej się pod stołem), oraz \( \int_0^1 B_2(x)B_3(x)dx = \int_0^1 B_1(x)B_2(x)dx \), tak naprawdę musimy więc policzyć jedynie cztery całki
\( \begin{bmatrix} \int_0^1 B_1(x)B_1(x)dx & \int_0^1 B_1(x)B_2(x)dx & \int_0^1 B_1(x)B_3(x)dx \\ == & \int_0^1 B_2(x)B_2(x)dx & == \\ == & == & == \\ \end{bmatrix} \)
Obliczamy jednowymiarowe całki z wielomianów
\( \int_0^1{B_1(x)B_1(x)dx} = \int_0^1 (\frac{1}{2} x^2)(\frac{1}{2} x^2)dx = \int_0^1 \frac{1}{4}x^4)dx = \frac{1}{4} (\frac{x^5}{5})|^1_0=\frac{1}{20} \)
\( \int_0^1{B_2(x)B_2(x)dx} = \int_0^1 (-x^2+x+\frac{1}{2})(-x^2+x+\frac{1}{2}) dx= \int_0^1 (x^4-2x^3+x+\frac{1}{4}) dx = \\ (\frac{x^5}{5})|^1_0-2(\frac{x^4}{4})^1_0+(\frac{x^2}{2})^1_0+(\frac{1}{4}x)|^1_0=\frac{1}{5}+\frac{1}{4}=\frac{9}{20} \)
\( \int_0^1{B_1(x)B_2(x)dx} = \int_0^1 (\frac{1}{2} x^2)(-x^2+x+\frac{1}{2})dx = \int_0^1 (-\frac{x^4}{2}+\frac{x^3}{2}+\frac{x^2}{4})dx = \\ (-\frac{x^5}{10})|^1_0+(\frac{x^4}{8})|^1_0 +(\frac{x^3}{12})|^1_0=-\frac{1}{10}+\frac{1}{8}+\frac{1}{12}=\frac{13}{120} \)
\( \int_0^1{B_1(x)B_3(x)dx} = \int_0^1 (\frac{1}{2} x^2)(\frac{1}{2}(1-x)^2)dx = \int_0^1 (\frac{x^4}{4}-\frac{x^3}{2}+\frac{x^2}{4})dx = \\ (\frac{x^5}{20})|^1_0 -(\frac{x^4}{8})|^1_0 +(\frac{x^3}{12})|^1_0=\frac{1}{20}-\frac{1}{8}+\frac{1}{12}=\frac{1}{120} \)
Uwzględniając symetrie obliczonych całek, wpisujemy je w stosowne miejsca w naszej matrycy i dostajemy
\( \begin{bmatrix} \int_0^1{B_1(x)B_1(x)dx} & \int_0^1{B_1(x)B_2(x)dx } & \int_0^1{B_1(x)B_3(x)dx} \\ \int_0^1{B_2(x)B_1(x)dx} & \int_0^1{B_2(x)B_2(x)dx } & \int_0^1{B_2(x)B_3(x)dx } \\ \int_0^1{B_3(x)B_1(x)dx} & \int_0^1{B_3(x)B_2(x)dx} & \int_0^1{B_3(x)B_3(x)dx } \\ \end{bmatrix} = \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} \\ \frac{13}{120} & \frac{9}{20} & \frac{13}{120} \\ \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \)
Sumując teraz wygenerowane w ten sposób kawałki macierzy dostaniemy dwie macierze tworzące nasz układ równań w postaci produktu Kroneckera macierzy.
Sumowanie przeprowadzamy w następujący sposób:
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{9}{20} +\color{red}{\frac{1}{20}}& \frac{13}{120} + \color{red}{\frac{13}{120}} & \color{red}{\frac{1}{120}} & \cdots \\ \frac{1}{120} & \frac{13}{120} +\color{red}{\frac{13}{120}} & \frac{1}{20} +\color{red}{\frac{9}{20}} +\color{blue}{\frac{1}{20}} & \color{red}{\frac{13}{120}} +\color{blue}{\frac{13}{120}} & \color{blue}{\frac{1}{120}} & \cdots \\ \cdots & \color{red}{\frac{1}{120}} & \color{red}{\frac{13}{120}} +\color{blue}{\frac{13}{120}} & \color{red}{\frac{1}{20}}+\color{blue}{\frac{9}{20}} +\color{green}{\frac{1}{20}} & \color{blue}{\frac{13}{120}} +\color{green}{\frac{13}{120}} & \color{green}{\frac{1}{120}} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ & \cdots & \color{green}{\frac{1}{120}} & \color{blue}{\frac{13}{120}} +\color{green}{\frac{13}{120}} & \color{red}{\frac{1}{20}}+\color{blue}{\frac{9}{20}} +\color{green}{\frac{1}{20} } & \color{red}{\frac{13}{120}}+\color{blue}{\frac{13}{120}} & \color{red}{\frac{1}{120}} \\ & & \cdots & \color{blue}{\frac{1}{120}} & \color{red}{\frac{13}{120}}+\color{blue}{\frac{13}{120}} & \frac{1}{20} +\color{red}{\frac{9}{20}} +\color{blue}{\frac{1}{20}} & \frac{13}{120} +\color{red}{\frac{13}{120}} & \frac{1}{120} \\ & & & \cdots & \color{red}{\frac{1}{120}} & \frac{13}{120} +\color{red}{\frac{13}{120}} & \frac{9}{20} +\color{red}{\frac{1}{20}} & \frac{13}{120} \\ & & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \)
Poza wypisanymi wyrazami macierzy umieszczonymi na pięciu przekątnych, pozostałe wyrazy macierzy są równe 0. Innymi słowy naszą matrycę przesuwamy wzdłuż przekątnej od lewego górnego rogu macierzy do prawego dolnego rogu, i sumujemy nakładające się na siebie wyrazy. W ten sposób uzyskamy dwie macierze, które tworzą za pomocą tak zwanego produktu Kroneckera naszą macierz ( 1 ).
Kolejny problem do rozwiązania to policzenie całek prawej strony. Całka prawej strony to próbkowanie bitmapy przemnożonej przez poszczególne funkcje testujące (B-spline'y używane do uśredniania bitmapy w obszarze na którym są określone).
Używamy od teraz nowej notacji, w której wprowadzamy dwuwymiarową funkcję B-spline, funkcję zmiennych \( (x,y) \), będącą iloczynem dwóch jednowymiarowych funkcji B-spline \( B_{k,l;2}(x,y)=B^x_{k;2}(x)B^y_{l;2}(y) \).
Zauważmy że przemnożenie bitmapy przez B-spline'a \( v(x,y)=B_{k,l;2}(x,y) \) oznacza iż całkę tą musimy policzyć jedynie po dziewięciu elementach na których określony jest B-spline \( B_{k,l;2}(x,y) \).
Musimy rozwiązać teraz dwa problemy. Pierwszy problem to jaki wzór ma dziewięć fragmentów naszego B-spline'a na tych dziewięciu elementach, na których jest on określony. Drugi problem, to jak policzyć całkę z bitmapy przemnożonej przez wielomian.
Zauważmy że jednowymiarowe B-spline'y składają się z trzech segmentów opisanych zgodnie z ( 1 ). W związku z tym, nasz dwuwymiarowy B-spline testujący na dziewięciu elementach na których jest określony składa się z odpowiednich iloczynów tych segmentów.
\( B_1(x)B_1(y), B_2(x)B_1(y), B_3(x)B_1(y) \\B_1(x)B_2(y), B_2(x)B_2(y), B_3(x)B_2(y) \\ B_1(x)B_3(y), B_2(x)B_3(y), B_3(x)B_3(y) \).
Jest to zilustrowane na Rys. 2.

Tak więc problem liczenia całki bitmapy przemnożonej przez funkcję testującą B-spline sprowadza się do problemu policzenia dziewięciu całek. Na przykład, jeśli chcemy przeliczyć całkę dla B-spline'a testującego \( B_{k,l;2}(x,y) \), wówczas musimy policzyć następujące dziewięć całek:
\( \int^1_0 \int^1_0 B_1(x)B_1(y) BITMAP((k-1+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_2(x)B_1(y) BITMAP((k+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_3(x)B_1(y) BITMAP((k+1+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_1(x)B_2(y) BITMAP((k-1+x)*maxx/N_x,(l+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_2(x)B_2(y) BITMAP((k+x)*maxx/N_x,(l+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_3(x)B_2(y) BITMAP((k+1+x)*maxx/N_x,(l+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_1(x)B_3(y)BITMAP((k-1+x)*maxx/N_x,(l+1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_2(x)B_3(y)BITMAP((k+x)*maxx/N_x,(l+1+y)*maxy/N_y) dxdy, \\ \int^1_0 \int^1_0 B_3(x)B_3(y)BITMAP((k+1+x)*maxx/N_x,(l+1+y)*maxy/N_y) dxdy. \)
W powyższych wzorach wyrażenia typu \( (k-1+x)*maxx/N_x \) oznaczają przeliczanie współrzędnych B-spline'a, którego całkujemy na indeksy pikseli bitmapy, na których on się rozpościera.
Zmienne x i y zmieniają się od 0 do 1. W obrębie jednego elementu, mamy \( maxx/N_x \) pikseli w kierunku \( x \)
oraz \( maxy/N_y \) pikseli w kierunku \( y \) , gdzie \( maxx,maxy \) oznacza rozmiar bitmapy (całkowitą liczbę pikseli w kierunku \( x \) i \( y \)) natomiast \( N_x,N_y \) oznacza liczbę elementów siatki w kierunku \( x \) i \( y \). Tak więc \( x*maxx/N_x \) zmienia się od 0 do \( maxx/N_x \), czyli obejmuje liczbę pikseli bitmapy w kierunku \( x \) na pojedyńczym elemencie, natomiast \( y*maxy/N_y \) zmienia się od 0 do liczby pikseli bitmapy w kierunku \( y \) na pojedynczym elemencie.
Dodatkowo przesuwamy nasz zakres tak żeby przeglądać piksele odpowiednich elementów, na których rozpięty jest nasz B-spline testujący. Na przykład wyrażenie \( (k-1+x)*maxx/N_x \) oznacza że przesuwamy się po pikselach elementu numer \( k \) w kierunku \( x \), \( (k+x)*maxx/N_x \) oznacza, że przesuwamy się po pikselach elementu numer \( k+1 \) w kierunku \( x \), \( (k+1+x)*maxx/N_x \) oznacza, że przesuwamy się po pikselach elementu numer \( k+2 \) w kierunku \( x \). Podobnie \( (l-1+y)*maxy/N_y \) oznacza, że przesuwamy się po pikselach elementu numer \( l \) w kierunku \( y \), \( (l+y)*maxy/N_y \) oznacza, że przesuwamy się po pikselach elementu numer \( l+1 \) w kierunku \( y \), oraz \( (l+1+y)*maxy/N_y \) oznacza, że przesuwamy się po pikselach elementu numer \( +2l \) w kierunku \( y \).
Jak więc policzyć takie całki po pikselach?
Musimy rozbić całkę na sumę całek po poszczególnych pikselach, na przykład dla pierwszej całki
\( \int^1_0 \int^1_0 B_1(x)B_1(y) BITMAP((k-1+x)*maxx/N_x,(l-1+y)*maxy/N_y) dxdy = \\ \color{red}{\sum_{i=1,maxx/N_x;j=1,maxy/N_y}}BITMAP((k-1)*maxx/N_x+i,(l-1)*maxy/N_y+j) \\ \color{blue}{\int_{(i-1,i)*\frac{1}{maxx/N_x}}} B_1(x)dx \color{blue}{ \int_{(j-1,j)*\frac{1}{maxy/N_y}}} B_1(y)dy \)
gdzie suma zaznaczona na czerwono rozpościera się po wszystkich pikselach z pojedynczego elementu (jest ich \( maxx/N_x * maxy/N_y \)), bitmapa jest próbkowana w punktach pikseli rozpostartych na danym elemencie, a całki zaznaczone na niebiesko, liczone są z funkcji B-spline po pojedynczych pikselach. Wystarczy teraz, że policzymy pojedyncze całki z trzech segmentów
\( \int B_1(x) dx = \int \frac{1}{2} x^2 dx = \frac{1}{2}\frac{x^3}{3}, \\ \int B_2(x) dx = \int (-x^2+x+\frac{1}{2}) dx = -\frac{x^3}{3}+\frac{x^2}{2}+\frac{1}{2}x, \\ \int B_3(x) dx = \int \frac{1}{2}(1-x)^2 dx = \int(1-2x+x^2)=x-x^2+\frac{x^3}{3} \)
i wstawimy je do wzorów na całki, na przykład
\( \color{blue}{ \int_{(i-1,i)*\frac{1}{maxx/N_x}} B_1(x)dx} = \left(\frac{1}{2}\frac{x^3}{3}\right)|^{ \frac{i}{maxx/N_x} }_{\frac{i-1}{maxx/N_x} } = \frac{1}{2} \left( \frac{ \left( \frac{i-1}{maxx/N_x} \right)^3- \left( \frac{i}{maxx/N_x}\right)^3 }{3} \right) \\ \color{blue}{\int_{(j-1,j)*\frac{1}{maxy/N_y}}B_1(y)dy } = \left( \frac{1}{2}\frac{x^3}{3} \right)|^{ \frac{j }{ maxy/N_y} }_{ \frac{j-1}{maxy/N_y} }= \frac{1}{2} \left( \frac{ \left( \frac{j-1} { maxy/N_y } \right)^3-\left(\frac{j}{maxy /N_y }\right)^3}{3}\right) \)
Liczenie całek funkcji testujących z bitmapy zrobiło się trochę skomplikowane. Podsumujmy, musimy policzyć wiele całek z bitmapy, przemnożone przez różne dwuwymiarowe B-spliny, tak zwane funkcje testujące. Każdą z tych całek liczymy tak że rozbijamy ją na całkę po dziewięciu elementach na których rozpięty jest B-spline testujący, i sumujemy te wynikowe dziewięć całek, wpisując otrzymaną liczbę do jednego wiersza wektora prawej strony. Z kolei żeby policzyć całkę po każdym z tych dziewięciu elementów na których rozpięty jest nasz B-spline testujący, całki te rozbijamy na poszczególne piksele które leża na tym elemencie, próbkujemy naszą bitmapę po tych pikselach i liczymy całkę z naszych segmentów B-spline'ów po wszystkich pikselach. Wszystko to razem sumujemy (suma po dziewięciu elementach na których leży B-spline testujący, suma po pikselach elementów).